The electronic structures, the equilibrium geometries and finite temperature 

properties of Na n (n=39-55) 
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Density-functional theory has been applied to investigate systematics of sodium clusters Na n in the 
size range of n= 39-55. A clear evolutionary trend in the growth of their ground-state geometries 
emerges. The clusters at the beginning of the series (n=39-43) are symmetric and have partial 
icosahedral (two-shell) structure. The growth then goes through a series of disordered clusters 
(n=44-52) where the icosahedral core is lost. However, for n >53 a three shell icosahedral structure 
emerges. This change in the nature of the geometry is abrupt. In addition, density-functional 
molecular dynamics has been used to calculate the specific heat curves for the representative sizes 
n= 43, 45, 48 and 52. These results along with already available thermodynamic calculations for n= 
40, 50, and 55 enable us to carry out a detailed comparison of the heat capacity curves with their 
respective geometries for the entire series. Our results clearly bring out strong correlation between 
the evolution of the geometries and the nature of the shape of the heat capacities. The results also 
firmly establish the size-sensitive nature of the heat capacities in sodium clusters. 

PACS numbers: 31.15.Ew, 31.15.Qg, 36.40.Ei, 36.40.Qv 



I. INTRODUCTION 

Physics and chemistry of clusters are very active ar- 
eas of research especially because of the emergence of 
nano science and nano technology* 1 Although major ef- 
forts have been spent into ground-state investigations, 
finite temperature properties are turning out to be very 
interesting. Such investigations are challenging, both ex- 
perimentally as well as theoretically. One of the first de- 
tailed measurements providing much impetus for theoret- 
ical work was on free sodium clusters by Haberland and 
co-workers.— These measurements reported the melting 
temperatures (T m ) of sodium clusters in the size range 
between 55 and 350 and remained unexplained for almost 
about a decade. 

The main puzzle was related to the irregular behav- 
ior of the melting temperature and the absence of any 
correlation between the peaks and the magic numbers ei- 
ther geometric or electronic. A good deal of simulation 
works has been carried out to explain the sodium data, 
most of the early work being with classical inter atomic 
potentials.—"^ It turned out that none of these could ob- 
tain qualitative and quantitative agreement with the ex- 
perimental data. Thus, it needed an ab initio density- 
functional method to achieve this. Indeed, much insight 
and excellent quantitative agreement has been obtained 
by density-functional molecular dynamics (DFMD) sim- 
ulations.£&LS&ia 

Recently a very different aspect of finite temperature 
behavior has been brought out by the experimental and 
later by theoretical work on gallium clustersi 11 ' 12 The 
experimental reports of Breaux and co- workers^ showed 
that in the size range of N—30 to 55, free clusters of gal- 
lium melt much above their bulk melting temperature. 
Interestingly, their experiment also showed that the na- 
ture of heat capacity is size sensitive. In fact, addition of 



even one atom changes the shape of the heat capacity dra- 
matically, e.g for Ga 30 and Ga 31 . A similar experimental 
observation has been reported for aluminum clusters in 
the size range n=49-62*i£ Such a size sensitive behavior 
has also been observed in DFMD simulation of Auig and 
Au2o^ A detailed analysis of the ground-state geome- 
tries of these clusters brought out the role of order and 
disorder in their geometries on the shape of the melting 
curve. A disordered system is shown to display a con- 
tinuous melting transition leading to a very broad heat 
capacity curve. However, the effect is subtle and descrip- 
tion of order and disorder needs careful qualifications. 

In spite of substantial experimental works on sodium 
clusters over a period of 10 years and or so, there is no 
firm and systematic evidence of size sensitivity. This is 
mainly due to the fact that the reported experimental 
data 2 - is for the size range of n=55-350 at discrete sizes. 
It is necessary to investigate the effect of addition of few 
atoms in a continuous manner in appropriate size range. 
However, it is not clear whether larger clusters having 
sizes of n >100 will also show this effect. 

An extensive ab initio study on the structural proper- 
ties of small Na clusters up to iV=20 has been reported 
by Rothlisberger and Andreoni The study reveals that 
pentagonal motifs dominate the structures above N=7. 
As expected most of the atoms in these clusters lie on 
the surface and a discernible core develops after about 
A/=15-16. The shapes after this sizes show signature of 
icosahedral structures. 

The finite temperature behavior of sodium clusters in 
the size range of n=8 to 55 has been reported; 7 The 
study reveals that it is not easy to discern any melt- 
ing peaks below n=25. However, the simulation data 
available at rather coarse sizes above n=40 already shows 
size-sensitive feature. In addition to this feature, our re- 
cent investigation^ on Nas7 and Na58 does bring out the 
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role of geometry and electronic structure on the melt- 
ing. Nevertheless, in order to draw definitive conclusions 
it is necessary to investigate the effect of addition of a 
few atoms in a continuous manner in appropriate size 
range. Therefore, we have chosen the size range of n=39- 
55 and have carried out detailed density-functional inves- 
tigations. The purpose of the present work is two fold. 
First, to obtain reliable equilibrium geometries for all the 
clusters in the size range of n=39-55 and to discern evo- 
lutionary trends. We note that Na4o has a symmetric 
partially icosahedron core and Nass is a complete icosa- 
hedron. Thus it is of considerable interest to examine 
the growth pattern from n=39 to n=55. The second 
purpose is to seek correlation between the nature of the 
ground-state and the evolutionary trends observed in the 
nature of their specific heats. Towards this end we have 
carried out extensive finite temperature simulations on 
representative clusters of size n= 43, 45, 48, and 52. 
Together with the already published results, this gives 
us access to specific heats for Na4o, Na43, Na 45 , Na48, 
Naso, Na52 and Na55, a reasonable representation across 
the series under investigation. Finally, we note that all 
the DFMD simulations reported so far have yielded ex- 
cellent agreement with the experimental data^^ These 
reports demonstrate the reliability of density-functional 
molecular dynamics in describing the finite temperature 
properties. 

The plan of the paper is as follows. In the next sec- 
tion (Sec. |n| we note the computational details. Sec. IIIII 
presents equilibrium geometries and their shape system- 
atic. Sec. IIVI presents the finite temperatures behavior of 
Na43, Na45, Na4 8 , Na 52 and finally we discuss the corre- 
lation between the ground states and nature of the spe- 
cific heats for all the available thermodynamics data. We 
conclude our discussion in Sec. [VJ 



simulation time of 2.5 ns per system. In order to get 
converged heat capacity curve especially in the region 
of co-existence, more temperatures were required with 
longer simulation times. We have discarded at least first 
30 ps for each temperature for thermalization. To ana- 
lyze the thermodynamic properties, we first calculate the 
ionic specific heat by using the Multiple Histogram (MH) 
technique. 23 ! 24 We extract the classical ionic density of 
states (Tt(E)) of the system, or equivalently the classi- 
cal ionic entropy, S(E) — fcglnO(i£), following the MH 
technique. With S(E) in hand, one can evaluate ther- 
modynamic averages in a variety of ensembles. We focus 
in this work on the ionic specific heat. In the canon- 
ical ensemble, the specific heat is defined as usual by 
C(T) = dU{T)/dT, where U(T) = J Ep(E, T) dE is 
the average total energy. The probability of observing 
an energy £ at a temperature T is given by the Gibbs 
distribution p(E, T) = n(E)exp(-E/k B T)/Z(T), with 
Z(T) the normalizing canonical partition function. Wc 
normalize the calculated canonical specific heat by the 
zero-temperature classical limit of the rotational plus vi- 
brational specific heat, i.e., Cq — (3N — 9/2)ks- 25 

We have calculated a number of thermodynamic indi- 
cators such as root-mean-square bond length fluctuations 
(<5rms), mean square displacements (MSD) and radial dis- 
tribution function (g(r)). The <5 rms is defined as 

A 2 ^ ((r 2 ,),-^) 2 ) 1 / 2 

drms " N(N - 1) 2. {rij)t ' W 

where N is the number of atoms in the system, is 
the distance between atoms i and j, and (. ■ -)t denotes 
a time average over the entire trajectory. MSDs for in- 
dividual atoms is another traditional parameter used for 
determining phase transition and is defined as, 



II. COMPUTATIONAL DETAILS 

We have carried out Born-Oppenheimer molecular dy- 
namics (BOMD) simulations 16 using Vanderbilt's ultra- 
soft pseudopotentials 17 within the local-density approxi- 
mation (LDA), as implemented in the VASP package. 18 
We have optimized about 300 geometries for each of 
the sodium clusters in the size range between n=39 and 
n=55. The initial configuration for the optimization of 
each cluster were obtained by carrying out a constant 
temperature dynamics simulation of 60 ps each at vari- 
ous temperatures between 300 to 400 K. For many of the 
geometries we have also employed basin hopping 19 and 
geneti o 20 ' 21 algorithms using Gupta potential^ for gen- 
crating initial guesses. Then we optimized these struc- 
tures by using the ab initio density-functional method^ 
For computing the heat capacities, the iso-kinetic BOMD 
calculations were carried out at 14 different temperatures 
for each cluster of Na4 3 , Na 45 , Na 48 , Na 52 in the range 
between 100 K and 460 K, each with the time dura- 
tion of 180 ps or more. Thus, it results in the total 



<*?(*)> = ]{f E t R ^om + 1) - R/(W)] 2 (2) 

TCI — 1 

where R/ is the position of the I th atom and we aver- 
age over M different time origins to m spanning over the 
entire trajectory. The interval between the consecutive 
to™ for the average was taken to be about 1.5 ps. The 
MSDs of a cluster indicate the displacement of atoms in 
the cluster as a function of time. The g(r) is defined 
as the average number of atoms within the region r and 
r + dr. 

We have also calculated the shape deformation param- 
eter (sdef), to analyze the shape of the ground state for 
all the clusters. The edef is defined as, 



£def = 




where Q x > Q y > Q z are the eigenvalues, in descend- 
ing order, of the quadrupole tensor 
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FIG. 1: The ground-state geometries of Na n (n=39-43). 
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FIG. 2: The ground-state geometries of Na n (n— 44-52). 



Qij — />'/, />'/ , 



(4) 



Here i and j run from 1 to 3, / runs over the number 
of ions, and Ru is the i th coordinate of ion / relative 
to the center of mass (COM) of the cluster. A spherical 
system (Q x — Q y = Q z ) has £d e /=l and larger values of 
EcZe/ indicates deviation of the shape of the cluster from 
sphericity. 
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FIG. 3: The ground-state geometries of Na n (n=53-55). 
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FIG. 4: The shape deformation parameter for Na n (n=39-55) 
as a function of cluster size. 



III. GEOMETRIES 

The lowest energy geometries of sodium clusters (Na„ 
n=39-55) are shown in Fig. Q] (n=39-43), Fig. [2] (n=44- 
52), and Fig. [3] (n=53-55). We have also plotted the 
shape deformation parameter Sdef and the eigenvalues 
of quadrupole tensor for the ground state geometries of 
these clusters in Fig. |4] and Fig. El respectively. 

It is convenient to divide these clusters into three 
groups. The clusters in the first group, shown in Fig. [T] 
are nearly spherical. The ground-state geometry for Na3g 
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FIG. 5: The eigenvalues of the quadrupole tensor for Na n 
(n=39-55) as a function of cluster size. 
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FIG. 6: The distance from the center of mass for each of the atoms ordered in the increasing fashion for Na n (n=39-55). The 
formation of the shells is evident from the sharp steps. 



is highly symmetric. This structure has three identical 
units, each based on the icosahedral motive and these 
three units are arranged as shown in the Fig. [TJ It is 
interesting to note that an addition of an extra atom 
changes the structure dramatically. It can be seen that 
the geometry of Na4o is based on icosahedral structure 
with missing 12 corner atoms and (111) facet as reported 
by Rytkonen et al £ An extra atom added to this struc- 
ture is accommodated near the surface and deforms the 
structure slightly. In addition to the deformation the dis- 
tance between the two shells is reduced by 0.3 A as com- 
pared to that in Na4o. A single atom added to Na4i is not 
accommodated in the structure, instead it caps the sur- 
face. However, the low-lying geometries for Na42 have a 
spherical shape without any cap (figure not shown). The 
lowest-energy structure of Na43 shows two caps symmet- 
rically placed on the opposite side of Na4i , accompanied 
by distortion of the icosahedral core. The second group in 
Fig- El consisting of the clusters with n=44-52 shows sub- 
stantial distortion of the icosahedral core and even loss 
of this core structure. These clusters essentially repre- 
sent the transition region from the two shell icosahedron, 
Na4o to three shell complete icosahedron, Na55. It can 
be seen that by adding one atom to Na43 the two shell 



core is destroyed. The growth from n=44 to n=52 shows 
successive stages of capping followed by rearrangement in 
the inner core. There is a dramatic change in the struc- 
ture as soon as we add one more atom to Na52- All the 
atoms rearrange to form an icosahedral structure as seen 
in Na53. Thus, the clusters in the last group namely, 
Na53 and Nas4 differ from a perfect icosahedron of Na55 
by an absence of two and one atom(s) , respectively (Fig. 

The nature of the changes in the shape of the clus- 
ters during the growth can be seen in Fig. 2] and Fig. 
[SJ The shape deformation parameter (sdef) increases to 
a value about 2 up to Na52 with slightly higher values 
for Na45 and Na4g (Fig. 2]). However, this value drops 
suddenly for Nas3. It is interesting to examine the three 
eigenvalues of quadrupole tensor (Q x , Q y , Q z ) shown in 
Fig. [5l It can be seen that two of the eigenvalues are 
nearly same up to Na52 while the third one continuously 
grows and indicates that the growth dominantly takes 
place along one of the directions. A prolate configura- 
tion has Q z S> Q x ~ Qy Thus, the majority of the 
clusters in the second group are prolate. The formation, 
the destruction and reformation of the shell structure is 
clearly seen in Fig. [S] In this figure we have plotted the 
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FIG. 7: (a) The normalized heat capacity and (b) the <5 rl 
for Na43. The peak in heat capacity curve is at 270 K. 
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FIG. 9: The radial distribution function calculated for Na43 
at five different temperatures. 



(a) (b) 

FIG. 8: Two low lying isomers of Na43. AE represents the 
energy difference with respect to the ground-state. 

distance of each atom from the center of mass arranged 
in the increasing fashion. Clearly small rearrangement of 
atoms yields a change in the structure from Na3g to Na4o- 
The two shell structures are observed till Na43. The for- 
mation of the shell structure is reflected in the formation 
of the sharp steps in the graph. As the size increases 
from Na43 to Na44 the shell structure is destroyed and 
seen again at Nas3. Thus, three shells begin to emerge 
at Na53. 



IV. THERMODYNAMICS 

We have calculated the ionic heat capacity and indi- 
cators like mean square displacements (MSD) and root- 
mean-square bond length fluctuations (<5 rms ) for four 
of the representative clusters in the investigated series 
which are Na43, Na45, Na 48 , and Nas2- 



We note that the thermodynamics of Na4(j£>£, Nas^i 
and Na5g££ has already been reported. Thus it is pos- 
sible to examine and analyze the systematic variations 
in the melting characteristic and correlate them with the 
equilibrium geometries across the entire range of sizes 
from 40 to 55. The heat capacity and <5 rms for Na43 are 
shown in Figs.[7{a) and[7jb), respectively. We also show 
typical low energy geometries (isomers) for Na43 m Fig- 
[5] The first isomer shown in Fig.[5Ja) has two atoms clos- 
est to each other capping the surface and second isomer 
shows a distorted shape and no caps. The heat capacity 
shows a weak peak around 160 K while the main peak 
occurs at 270 K. An examination of the motion of ionic 
trajectory seen as a movie indicates that isomcrization 
(Fig. Ufa)) is responsible for the weak peak. 

It is interesting to observe the changes of radial dis- 
tribution function (RDF) as a function of temperature 
which is shown in Fig. [5] At low temperatures the shell 
structure is clearly evident. The pattern seen at 180 K 
and 210 K are mainly due to the fluctuation of the cluster 
between the ground-state and low-lying states. At 300 K 
and above the RDF shows the typical melting behavior 
of a cluster. The <5 rms in Fig.[7](b) shows the effect of iso- 
merization around 160 K. It can be seen that the melting 
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FIG. 10: (a) The normalized heat capacity and (b) the <5 rms 
for Na45 . The heat capacity curve shows a two stages melting 
process. 
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FIG. 12: (a) The normalized heat capacity and (b) the <5 rl 
for Na4g, (c) The normalized heat capacity and (d) the <5 rl 
for Na 5 2. 



FIG. 11: The MSDs of individual atoms calculated for Na45 
(a) at 210 K and (b) at 225 K over the last 150 ps. 



region is of the order of 60 K. 

Let us recall that Na45, Na48 and Na52 belong to the 
second class of disordered clusters. Among these Na45 
shows some signature of partial shells. The heat capacity 
and 5 rms for Na45 are shown in Figs. [TUTal and [TUTb). 
respectively. The heat capacity of Na45 shows a first 



peak around 230 K and a second peak around 300 K. We 
also show the MSDs of individual atoms at 210 and 225 K 
in Figs. [TIT a) and fTTT b). respectively. We have observed 
the ionic trajectories as a movie in the temperature at 
225 K. It turns out that about one third of atoms in the 
cluster "melt" at this temperature. This is brought out 
by the contrasting behavior of the MSDs as shown in 
Figs. [TTTa) and [TTT bV Interestingly, all these 15 atoms 
are on the surface, indicating that surface melting takes 
place first. The peak around 225 K in the heat capacity 
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FIG. 13: The normalized heat capacity as a lunction of tem- 
perature for Nan, n=40, 43, 45, 48, 50, 52, and 55. 



is due to partial melting of these surface atoms. We show 
the <5 rms as a function of the temperature in Fig. [Tfflb) . 
There is a sharp increases in <5 rms around 210 K and a 
slow rise after 240 K consistent with the behavior of the 
heat capacity. Thus in Na45 melting takes place in two 
stages over the range of 120 K. 

The heat capacity and <5 rms for Na48 an d Na52 are 
shown in Fig. [T3J The heat capacity of Na4§ (Fig. [T27a)) 
is very broad indicating almost continuous phase change 
starting around 150 k. Thus it is difficult to identify a 
definite melting temperature. A similar behavior is also 
seen in that of Nas2 except for the small peak seen 180 K 
due to isomerization (Fig.[T27c)). It is interesting to note 
that Srms for both clusters also show a gradual rise from 
about 150 K to 350 K indicating a continuous melting 
transition as shown in Figs. [T2T b) and [T2T d). 

The systematic evolution of heat capacities can be bet- 
ter appreciated by examination of all the available data 
(calculated from density- functional simulation). In Fig. 
[TSlwe show the specific heats for all the available clusters 
between Na4o and Nass. The most symmetric cluster 
Na55 shows the sharpest peak in the heat capacity. The 
heat capacity of Na4o aud Na43 (partial icosahedral struc- 
tures) have well recognizable peaks which are broader 
than that of Nas5. The disordered phase of the growth 
is clearly reflected in the very broad heat capacities seen 
around n=50. 



V. SUMMARY AND CONCLUSION 

The ab initio density-functional method has been ap- 
plied to investigate systematic evolutionary trends in the 
ground state geometries of the sodium clusters in the size 
range of n=39-55. The DFMD finite temperature simula- 
tions have been carried out for representative clusters. A 
detailed comparison between the heat capacities and the 
geometries firmly establishes a direct influence of the ge- 
ometries on the shapes of the heat capacity curves. The 
heat capacities show size sensitivity. The growth pattern 
shows a transition from ordered — > disordered — > or- 
dered sequence. The corresponding heat capacities show 
a transition from peaked to a very broad to peaked se- 
quence. It is seen that addition of a few atoms changes 
the shape of heat capacity very significantly. We believe 
that the size sensitive feature seen our simulation is uni- 
versal. It may be noted that such a feature has been 
observed experimentally in the case of gallium and alu- 
minum and in the case of DFMD simulations for gold. 
We await the experimental measurements of the heat ca- 
pacities on the sodium clusters in these range showing 
the size sensitivity. 
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